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Abstract 

The alignment of Saturns magnetic pole with its rotation axis precludes the use of magnetic 
field measurements to determine its rotation period (ll. The period was previously determined 
from radio measurements by the Voyager spacecraft to be lOh 39m 22.4s ||2l. When the Cassini 
spacecraft measured a period of lOh 47m 6s, which was additionally found to change between 
sequential measurements SIH |5l, it became clear that the radio period could not be used to 
determine the bulk planetary rotation period. Estimates based upon Saturn’s measured wind 
fields have increased fhe uncerfainfy even more, giving numbers smaller fhan fhe Voyager ro- 
fafion period, and af presenf Safum’s rofafion period is fhoughf fo be befween lOh 32m and 
lOh 47m, which is unsalisfaclory for such a fundamenlal properly. Here we report a period 
of lOh 32m 45s ± 46s, based upon an oplimizalion approach using Salurns measured gravila- 
fional field and limils on fhe observed shape and possible infernal density profiles. Moreover, 
even when solely using fhe consfrainfs from ifs gravifafional field, fhe rofafion period can be 
inferred wilh a precision of several minutes. To validafe our mefhod, we applied fhe same 
procedure fo Jupifer and correcfly recovered ifs well-known rofafion period. 
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Previous theoretical attempts to infer Satums rotation period have relied on wind observations 
derived from cloud tracking at the observed cloud level iH . One theoretical approach was based 
on minimizing the 100 mbar dynamical heights [|71 with respect to Saturns measured shapeS, while 
a second approach was based on analyzing the potential vorticity in Saturns atmosphere from its 
measured wind profile [(81. The derived rotation periods were found to be lOh 32m 35s ±13s, 
and lOh 34m 13s ± 20s, respectively. Our optimization method is based on linking the rotation 
period of Saturn with its observed physical properties and their uncertainties, in particular, the 
gravitational field. The method allows us to derive Satums rotation period for different types of 
constraints, and does not rely on a specific interior model, equation of state, wind properties, or 
other indirect measurements. 

Previous theoretical investigations to infer Saturn’s rotation period have relied on wind ob¬ 
servations derived from cloud tracking at the observed cloud level 0 . One theoretical approach 
was based on minimizing the 100 mbar dynamical heights flOl with respect to Saturn’s measured 
shape [|71, while a second approach was based on analyzing the potential vorticity in Saturn’s 
atmosphere from its measured wind profile [[ 8 l. The derived rotation periods were found to be 
lOh 32m 35s ±13s, and lOh 34m 13s ± 20s, respectively. 

The gravitational moments and the internal density profile can be related through the small¬ 
ness parameter m = oj'^R^/GM, where R is the planet’s mean radius, M is its mass, G is 
the gravitational constant, and oj = 27r/P is the angular velocity associated with the rotation 
period P [fTOl [TTI . The even gravitational moments can be expanded as a function of m by 
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J 2 n = X] whcrc a 2 n,k STC cocfficients that are determined by the radial density distri- 

k=n 

bution (see Methods). The expansion can go to any order of n; since at present only J 2 , J 4 and Jq 
are known for Saturn (and Jupiter), in this study we take n = 1,2,3. 

The relation for J 2 n shows that the measured gravitational moments are determined from the 
combination of the internal density distribution {a 2 n,k) (^nd rotation (m). Our goal is to find a solu¬ 
tion for m and a 2 n,k that minimizes the difference between the observed and calculated J 2 n within 
the observed uncertainties. The a 2 n,k can be expressed by a combination of figure functions (see 
Methods) that represent a given internal density profile, and can then be linked to the gravitational 
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harmonics [fTOl W2\. [T3]| . However, sinee in this case there are only three equations and seven un¬ 
knowns (six figure funetions and the smallness parameter), there is no unique solution. As a result, 
the solution is found by using a statistieal optimization approaeh. 

We define an optimization funetion as the sum of the normalized absolute differenees between 
the observed gravitational moments and the ealeulated gravitational moments, given by 



where J 2 n are the calculated moments, are the measured moments, and are the measure¬ 
ment uneertainties of the measured gravitational moments [fT4l [T5ll . The optimization proeedure 
begins with an initial guess of the various parameters being randomly spread throughout the physi- 
eal bounds of eaeh parameter. This is repeated 2000 times to aehieve statistieal signifieanee. From 
these 2000 eases we eompute the rotation period and its standard deviation (see Methods). An 
example of the derived solutions using our optimization method is presented in Fig. 1. 

The entire set of solutions for Saturn are summarized in Fig. 2 whieh shows Peak (dots) and 
its Icr standard deviation (blue-shading). We first present solutions that are completely uneon- 
strained in radius and density strueture, and where the rotation period (gray-shading) is allowed 
to vary widely. The faet that the ealeulated standard deviation is mueh smaller than the allowed 
range (blue shading being mueh narrower than the gray shading) indieates that knowledge of the 
gravitational moments ean be used to narrow the possible range of rotation periods. In addition, 
as the initial range of the possible rotation periods is narrower, the derived rotation period ean be 
determined with higher accuracy. For the smallest range in rotation period (left dot in Fig. 2a) we 
derive a rotation period of lOh 43m 10s ± 4m. The fact that the uneertainty in rotation period is 
deereased significantly without enforcing tight constraints on the model emphasizes the strength 
of this method. Nonetheless, without any eonstraints on the shape the solution for the rotation pe¬ 
riod still has a relatively large range of solutions. In reality, oeeultation measurements ll^ fT^fTTll 
provide bounds on the shape of the planet (radius vs. latitude), and as shown below this allows to 
further eonstrain the rotation period [|7l [T8l[T9l . 

The best measurement uncertainty of Saturn’s radii from radio and stellar oeeultation is ~6 km 
liTTIl although the actual uncertainty eould be larger due to the unknown eontribution of the atmo- 
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spheric dynamics to the measured shape [|20l . We therefore explore a range of uncertainty in mean 
radius between 6 and 80 km. The results for this case are shown in Fig. 2b where Pcaic and its 
standard deviation vs. the uncertainty in observed radius Pobs are shown. The standard deviation 
(blue-shading) of Pcaic decreases with decreasing uncertainty in the radius. For an uncertainty of 
6 km in Saturn’s mean radius we derive a rotation period of lOh 34m 22s ± 3.5m. It is clear that 
the parameter space of possible solutions narrows when the constraint of Saturn’s measured mean 
radius is included. Yet, geopotential variations due to atmospheric dynamics affect the shape of 
the planet, and therefore caution should be taken when considering these measurements. By taking 
this hierarchal approach we are able to isolate the uncertainty given estimates of shape and internal 
structure separately. More conservative uncertainties in radius (tens of kilometers) yield longer 
rotation periods, thus giving solutions closer to the Voyager rotation period (see Fig. 2b). 

The uncertainty in Peak can be decreased even further if we also limit the range of the figure 
functions, i.e., the density profile (Fig. 2c). Limiting the figure functions to within a range implied 
by interior structure models (see Methods), the derived period is found to be lOh 32m 45s ± 46s. 
This rotation period is in agreement with previous calculations that derived Saturn’s rotation pe¬ 
riod by using a fit to its measured shape E HHl. The fact that the rotation period is shorter than 
the Voyager rotation period also implies that the latitudinal wind structure is more symmetric, thus 
containing both easterly and westerly jets as on Jupiter ||8l. Although the smallest possible uncer¬ 
tainty in rotation period is desirable, there is a clear advantage in not specifying constraints on the 
density profile, and keeping the method as general as possible. 

Different from Saturn, the rotation period of Jupiter is well determined due to its tilted mag¬ 
netic field. Jupiter’s measured rotation period (system III) is 9h 55m 29.69s [|2T1 l22l . In order 
to verify the robustness of our results we apply this method also for Jupiter (Fig. 3). When only 
the gravitational moments are used as constraints (Fig. 3a), like for Saturn, the uncertainty in the 
calculated rotation period is much smaller than the allowed range, and converges toward Jupiter’s 
rotation period. Fig. 3b shows the sensitivity of the derived period, when the uncertainty in period 
is ± 0.5 hr around the measured value, for a range of possible mean radii. Like for Saturn, the stan¬ 
dard deviation of Peak decreases with decreasing AR. When the variation in Pobs is taken to be 
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6 km, a rotation period of 9h 56m 6s ± 1.5m is derived, eonsistent with Jupiter’s measured rotation 
period. When we also add eonstraints on the figure funetions, the derived rotation period beeomes 
9h 55m 57s ± 40s, showing that our method reproduees Jupiter’s rotation period suoeessfully. 

The determination of Saturn’s J 2 n is expeeted to improve signifieantly following Cassini’s end- 
of-mission proximal orbits. To test whether a more aeeurate determination of Saturn’s gravitational 
field will allow to better eonstrain its rotation period, we repeat the optimization with the expeeted 
new uncertainty on the gravitational moments (A J 2 n ~ 10“®) [|^ around the currently measured 
values. The solution making no assumptions on the density profile is shown in Fig. 4a. Since 
Jupiter’s gravitational field will be more tightly determined by Juno [|^ l24l we do a similar anal¬ 
ysis for Jupiter (Fig. 4b). While for Jupiter the calculated rotation period remains the same with 
the more accurate gravitational field, for Saturn the calculated uncertainty of Peak decreases by 
~15%. We therefore conclude that the future measurements by Cassini will be important for fur¬ 
ther constraining Saturn’s rotation period. 
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Smallness parameter m 


Figure 1: An example of the statistical distribution of solutions for Saturn’s rotation period. 

For this specific case, the initial possible range of rotation periods is taken to have 0.5 h uncer¬ 
tainty around the Voyager radio period (hereafter, Pvoy)- The calculated mean radius was set to be 
within 20 km of Saturn’s observed mean radius. The solution is based on an ensemble of 2000 in¬ 
dividual sub-cases, each of them representing a case with specific random initial conditions within 
the defined parameter space, a. A scatter plot of the distribution of solutions on the plane of 
the calculated rotation period Peak minus Pvoy = lOh 39m 22s and AJ 2 . Each blue dot repre¬ 
sents one sub-case solution. The inner and outer black circles show the first and second standard 
deviations, respectively, b. The distribution of the derived rotation period as a function of small¬ 
ness parameter m and the calculated mean radius Peak minus the observed mean radius of Saturn 
(Pobs = 58, 232 km). 




































Saturn 



J AP [hr] 



c 



Figure 2: Solutions for Saturn’s rotation period, a. The calculated period Peak (blue dots) and 
its Icr standard deviation (blue shading) for a large range of cases for which the assumed possible 
range in rotation period varies between 0.25 h and 5.5 h (gray shading) around P^oy (black-dashed 
line), b. Peak and its la standard deviation (blue shading) using AP = 0.5 h vs. the assumed 
uncertainty in Saturn’s observed mean radius Pobs- c. Same as (b) but when the figure functions 
are also constrained. 
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Figure 3: Solutions for Jupiter’s rotation period, a. The ealeulated period Pcaic (blue dots) and 
its Icr standard deviation (blue shading) for a large range of oases for whioh the assumed possible 
range in rotation period varies between 0.25 h and 5.5 h, i.e., between ~ 5 h and 15 h (gray shading) 
around Jupiter’s measured period (black-dashed line), b. Peak and its la standard deviation (blue 
shading) using AP = 0.5 h vs. the uncertainty in the assumed uncertainty in Jupiter’s observed 
mean radius Pobs- c. Same as (b) but when the figure functions are also constrained. 
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Figure 4: Solutions for the rotation periods of Saturn (a) and Jupiter (b) when assuming 
improved gravity data. Shown are Pcaic and its Icr standard deviation (blue shading) when 
setting assuming A J 2 n = 10“® and AP = 0.5 h. The calculated period is given vs. the assumed 


uncertainty in the observed mean radius Pobs- 
















1 Methods 


1.1 The theory of figures 

The theory of figures was first introdueed by Clairaut [[25ll who derived an integro-differential 
equation for calculating the oblateness of a rotating planet in hydrostatic equilibrium with a non- 
uniform density profile. The method was further developed by Zharkov & Trubitsyn [flOl who 
presented a theoretical description to connect the density profile of a hydrostatic planet with its 
gravitational moments J 2 n, extending the theory to an arbitrary order. The basic idea of the method 
is that the density profile of a rotating planet in hydrostatic equilibrium can be derived by defining 
the layers as level surfaces., i.e., surfaces of a constant potential (called the effective potential) that 
is set to be the sum of the gravitational potential and the centrifugal potential [fTOl fTll 

u = (“) J 2 nP 2 n (cos 6*)^ + sin^ 6*, (2) 

where r is the radial distance, a is the equatorial radius of the geoid, GM is its mass multiplied by 
the gravitational constant, 0 is the colatitude, and oj is the angular velocity given by 27r / P, with 
P being the rotation period. The internal density profile and the gravitational moments are linked 
through the smallness parameter m = oo'^R^/GM, where R is the mean radius of the planet. The 
gravitational moments can be expanded as a function of m by, 


J 2 = 77102,1 + 771^02,2 + 771^02,3 

(3a) 

J4 = 777^04,2 + 777^04,3 

(3b) 

Je = 777^06,3 

(3c) 


where a 2 n,k are the expansion coefficients in smallness parameter. As J 2 S> J 4 S> Jq higher 
order harmonics correspond to a higher order expansion in m. The gravitational coefficients J 2 n 
are determined from the combination of the internal density distribution and rotation period. As a 
result, unless the density profile of Saturn (or any other giant planet) is perfectly known there is no 
simple way to derive the rotation period and vice versa. 




For the investigation of planetary figures, the equation for level surfaees ean be written in 
the form of a spheroid that is a generalized rotating ellipsoid. Then, the planetary radius r at 
every latitude ean be expressed as a funetion of the polar angle d (eolatitude), and the flattening 
parameters /, k and h by ifTOlfT^ : 


r {9) = a 


1 — / cos^ 9 — sin^ 29 + ^ ^ 


sin^ 29 


(4) 


where f = {a — b)/a is the flattening (with b being the polar radius), and k and h are the seeond- 
order and third-order eorreetions, respectively [fTOll . While / is strictly the flattening of the object, 
k and h represent the departure of the level-surface from a precise rotating ellipsoid to second and 
third order in smallness parameter, and their values are expected to be much smaller than /. 

To third order, the three flattening parameters /, k, h at the planetary surface (the effective 
potential surface) can be written as sum of figure functions defined by ifT^ . 


/ = mFi -I- m^F 2 -f (5a) 

k = m^K2 + (5b) 

h = m^H^. (5c) 

Finally, using the relation between the first three even gravitational coefficients and the figure 
functions for a density profile that is represented by a 6 th order polynomial [|7l [l 2 l [T3l and by 
applying the theory of figures as a set of differential equations, the gravitational coefficients and 
the figure functions can be related as power series in the small rotational parameter m (see equation 
72 in Schubert et al. 2011 ifT^ l. Since only J 2 , J 4 and Jq are currently known for Saturn (and 
Jupiter) we expand only up to third order in m. Although higher-order harmonics are not expected 
to be zero, the corrections will be 0 (m"^) and therefore their contribution will be small. 


1.2 The optimization method 

Since the flattening parameters (and figure functions) depend on the density distribution that is 
unknown, we take a general approach that is designed to relate the planetary rotation period to 





its gravitational field without putting tight constraints on the internal structure. We therefore de¬ 
veloped an optimization method that searches for the solutions that reproduce Saturn’s measured 
gravitational field within the widest possible pre-defined parameter space. The figure functions 
(F, if, K) are allowed to vary over their widest possible physical range, and the smallness pa¬ 
rameter m is allowed to vary within a range that reflects the uncertainty in the rotation period P. 
A solution for these parameters is sought while meeting the requirement that Saturn’s measured 
physical properties are reproduced. 

First, an optimization function is defined as the sum of the normalized absolute differences 
between the observed moments and the calculated moments and is given by: 



where J 2 , J 4 and Jg are the gravitational coefficients calculated using Eqs. 2a-2c, and 

are the measured gravitational moments, and Aand AJg’’® are the uncertainties 
on the measured gravitational coefficients lfT4l [TSll . Since the observations include only the first 
three even harmonics everything is calculated to third order, but the method can be modified to 
include higher order terms. The data that are used by the model are summarized in Extended Data 
Table 1. 

Next, we minimize the optimization function Y with respect to the control variables Fi, F 2 , F 3 , 
K 2 ,K^, and m, i.e., the figure functions, and the smallness parameter. Starting from an ar¬ 
bitrary initial guess for each of the seven control variables (within the predefined limits), a so¬ 
lution is sought such that the optimization function reaches a minimum. Several nonlinear con¬ 
straints are imposed while searching for the solution. Eirst, we require that the difference be¬ 
tween each calculated and the measured gravitational coefficients must be smaller than the uncer¬ 
tainty of the measurement error, i.e., | J 2 — < 0 , | J 4 — — |AJ 4 *’®| < 0 , and 

\Jq — 0 ^ 6 I — I< 0. Note that this requirement is additional to the minimization of Y since 
we ask that not only the overall difference between the observed and calculated gravity moments 
is minimized, but that individually, each of the calculated moments stays within the uncertainty of 
its observed counterpart. 

The parameter / is the planetary flattening, as a result, / must be a small positive number (for 





Saturn / ~ 0.1). The second and third order corrections, k and h, are significantly smaller than 
/, but could be both positive and negative. Thus, in order to keep our calculation as general as 
possible we let the three flattening parameters to vary between their maximum physical values, -1 
and 1. In Extended Data Fig. 1 we show the calculated values for /, h, k for Saturn for the case in 
which the figure functions are not constrained and Ai? is taken to be 50 km. / is found to be of 
the order of 0.1, consistent with the measured flattening of Saturn OEQI, while the second-order 
and third-order corrections are found to be of the order of 10“^. As Fi is the first order expansion 
for /, / — Fim = 0{m^), meaning that for Saturn |Fi| < 1. Similarly expanding recursively the 
other coefficients of /, and also for k and h, implies that all figure functions are bound between -1 
and 1. In order to keep our calculation as general as possible we let all the figure functions to vary 
between -1 and 1. The solution though, which must also fit the gravitational field, constrains the 
flattening parameters and figure functions to a much narrower range. The solution is derived by 
using a numerical algorithm that is designed to solve constrained nonlinear multivariable functions. 
We use a sequential quadratic method that formulates the above nonlinear constraints as Lagrange 
multipliers Ii26l . The optimization is completed once the tolerance values for the function (10“^) 
and the constraints (10“^^) are met. 

A single optimization would be sufficient if the problem was well defined. In such a case, there 
would have been a unique solution that is independent of the initial guess. However, since in our 
case there are only three equations (Eqs. 2a-2c) and seven unknowns (six figure functions and the 
smallness parameter), the problem is inherently ill-defined and therefore has no unique solution. 
Nevertheless, we can still reach a solution using a ’’statistical” approach in which we repeat the 
optimization process enough times to achieve a statistically stable solution. In each case, the initial 
guess of the various parameters is chosen randomly within the defined bounds of each parameter. 
A statistical significance is reached when we repeat the optimization 2000 times (verified with 10^ 
long optimizations). We can then use the solutions from the 2000 optimizations to compute the 
mean value and its standard deviation for each variable. An example for a specific case is given 
in Fig. 2 of the Extended Data where the solutions for the gravitational moments (Fig. 2a-c in 
Extended Data) are distributed around the mean value, and the distribution of the solutions for the 


figure functions (Fig. 2d-j in Extended Data) has a large range. From each such an experiment we 
eventually calculate two numbers: the mean rotation period Pcaic and its standard deviation. 

1.3 Saturn’s rotation period with an improved measurement of its gravita¬ 
tional field 

Another objective of this research is to understand whether the improved gravity measurements 
of the low-order gravitational moments (J 2 , J 4 , Jg) by Cassini Solstice mission can be used to 
better constrain Saturn’s rotation period. Since the data are not yet available, but are expected to 
be within the current uncertainty, all we can do it to estimate the rotation period and its standard 
deviation when the uncertainty on the gravitational harmonics (A J 2 n) is of the order of 10 “® while 
using the currently known values of J2, J4, and Jg. The result for this exercise is presented in 
Fig. 4 of the main paper. It is found that for Saturn this yields a 15% improvement in the derived 
standard deviation of its rotation period. However, it is important to remember that the true values 
of the gravitational moments can be any value within the current uncertainty. We find that the 
order of magnitude of the standard deviation is not very sensitive to the actual value of J 2 n but is 
more affected by the allowed uncertainty (i.e., A J 2 n), we can therefore conclude that the Cassini 
proximal orbit measurement is useful for further constraining Saturn’s rotation period. 

1.4 Accounting for the planetary shape 

Our optimization method can include additional constraints. Since the planetary shape could be 
used to constrain the rotation period [fT^| 20 l , we also run cases in which we account for Saturn’s 
shape (see Eq. 4). Then the optimization includes the constraint that the calculated mean radius R 
should be consistent with the mean radius that is inferred from measurements. Thus, the calculated 
mean radius Peak should be less than a specified uncertainty, i.e., 

|Pcalc - Pobsl - |AP| < 0, (7) 

where Pobs is the mean radius estimated from measurements of the planetary shape, and AP is the 
uncertainty associated with the measured radius. In the standard case we set this uncertainty to be 


40 km, which is large compared to the measured uncertainty in Saturn’s shape [|9]|. This provides 
a fourth equation to our optimization method and allows a considerable reduction to the rotation 
period uncertainty (Fig. 2b and Fig. 3b in the main paper). 

Although Saturn’s measured shape (radius as a function of latitude) is well determined from 
occultation measurements ||9l[T^, one should note that there is a difference between the measure¬ 
ment uncertainty, that is estimated to be ~6 km [T3 and the actual uncertainty that is of the 
order of a few tens of kilometers [l20ll . The actual uncertainty is relatively large because the plan¬ 
etary measured shape is also affected by atmospheric winds which distort the hydrostatic shape. 
The equatorial region of Saturn is affected by the large equatorial winds, and indeed the dynam¬ 
ical heights of the equator are found to be 120 km dH[IS|7l|^. On the other hand, the polar 
region is less affected by winds, and therefore the polar radii better reflects Saturn’s hydrostatic 
shape. There are, however, no available occultation measurements of Saturn’s polar regions. In 
addition, Saturn’s north-south asymmetry in wind structure introduces an additional uncertainty 
in determining its polar radius dll. As a result, a more conservative uncertainty in Saturn’s mean 
radius is estimated to be ~40 km d20l . 

Interestingly, the solution for the rotation period for the case without the shape constraint does 
not necessarily contain the solution when the constraint on the shape is included. This is caused by 
the fact that taking into account only the gravitational moments, for Saturn, leads to a solution with 
relatively long rotation periods, while the measured shape pushes to shorter rotation periods. This 
effect is illustrated in Extended Data Fig. 3 where the solutions for the rotation period is shown in 
the phase-space of the constraints for AP and AR. When there is no constraint on the shape, and 
AP is large (upper-right, red-region), the solution converges to a relatively long rotation period; 
as AP decreases, solutions with shorter rotation periods can be found (upper-left, blue region, see 
also figure 2a in main text). When the constraint on the shape is included, even when the range of 
the rotation period is large (blue region, bottom right), the solution converges into a short rotation 
period. The dashed line shows the transition between the region where the constraint on the period 
(above the dashed line) to the regime where the constraint on the shape is more important. For the 
relevant physical region we are interested in, as the constraint on the shape is increased (Ai?obs 


decreases), it becomes more important than the constraint on the rotation period and therefore it 
becomes more dominant leading to a shorter rotation period and outside of the range of the solution 
without the shape constraint. This behavior does not exist in the case of Jupiter as can be seen from 
Figure 3a in the main text. For Jupiter, the solutions reproduce the measured rotation period even 
for large ranges in AP and AR, and the two constraints are consistent with the same rotation 
period. 

1.5 Constraining the figure functions 

We also present results for cases in which the figure functions are constrained as well (Fig. 2c and 
Fig. 3c in the main paper for Saturn and Jupiter, respectively). In these cases, the figure functions 
are limited to be within a range that is determined from realistic interior models. In order to put 
limits on their values, we run two limiting interior models for both Saturn and Jupiter and derive 
the values of the figure functions. The first case is one of a massive core for which we assume a 
constant-density core with a density of ~1.5xl0^ kg m“^ reaching 20% of the planets radius. In 
the second case, the density is continuous and is represented by a 6th order polynomial. For this 
case, the first degree term of the polynomial is missing so that the derivative of the density goes 
to zero at the center. Another constraint sets this value to zero at the core-envelope boundary for 
models with cores. We then use the derived values of the figure functions for each case to limit the 
values of the variables Fi, F 2 , F 3 , K 2 , H^. 

The figure functions we derive for the massive core and continuous density profile cases for 
Saturn and Jupiter are summarized in Extended Data Table 2. The two density profiles for Saturn 
(top) and Jupiter (bottom) are shown in Extended Data Eig. 4. The values of the figure functions for 
interior models intermediate to these extreme cases all lie within the values of the extreme cases, 
implying that we have taken an exclusive set of values. In order to make sure that we account for 
a relatively large range of possible interior models, even within this fairly constrained case, we 
allow the figure function values to vary around the average value between the two models by a 
factor of two of the difference. The results (for Saturn) are shown in Extended Data Eig. 5. It is 
clear that in this case, due to the limitation on the figure functions, the parameter space of possible 



solutions is smaller allowing a more aeeurate determination of Pcaic- This range still aeeounts for 
a large variation in Saturn’s density profile. While there is an improvement in the determination 
of Pcaic when the shape and figure funetions are tightly eonstrained, there is also a elear advantage 
in keeping the method as general as possible. The inferred result is then not assoeiated with a 
speeifie interior model and/or does not relay on shape measurements. Our method is therefore also 
useful for estimating the rotation period of giant planets with less aeeurate determinations of their 
physieal properties. For the iey planets, Uranus and Neptune, only J 2 and J 4 are eurrently known, 
yet a simplified version of this optimization (to seeond order) ean be applied and gives rotation 
periods within 2% of the Voyager radio periods, allowing an independent method for estimating 
their rotation periods from their gravitational moments. Furthermore, our method eould also be 
applied for deriving the rotation periods of exoplanets for whieh the gravitational moments ean be 
estimated [l27l|28l . 
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Figure 5: Extended Data. The calculated flattening parameters by the model without con¬ 
straining the figure functions. Shown are /, fc, h for Saturn with Ai? = 50 km. 








Figure 6 : Extended Data. An example from our statistical optimization model for deriving 
the rotation period. The results are shown for a case for which the range of rotation period is 
lOh 24m - lOh 54m. The solution is based on a combination of 2000 individual sub-cases, each 
of them representing a case with specific random initial conditions within the defined parameter 
space, a-c A scatter plot (similar to Fig. la) of the distribution of solutions on the plane of the 
calculated rotation period Pcaic niinus Pvoy and A J 2 , A J 4 , A Jg, respectively. Each blue dot rep¬ 
resents one sub-case converged solution. In a the inner and outer black circles show the first and 
second standard deviations, respectively, d-j - the distribution of solutions for the figure functions 
A2, A's, P3, Pi, P2, P3, respectively. 
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Figure 7: Extended Data. Saturn’s calculated rotation period vs. the uncertainty in the as¬ 
sumed rotation period and radius. Shown Pcaic as a function of AP [hr] and AR [km]. The 
dashed line presents the transition between the region where the constraint on the rotation period 
(above the dashed line) to the regime where the constraint on the shape is more dominant. 
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Figure 8: Extended Data. Radial density profiles for two different interior models for Sat¬ 
urn (top) and Jupiter (bottom). The black curves correspond to models with very large cores 
while the blue curves are no-core models in which the density profile is represented by 6th order 
polynomials. For the massive-core case we assume a constant core density of ~ 1.5x10'^ kg m ^ 
reaching 20% of the planets radius. The density profiles are constrained to match the planetary 
mass, J2, J4, J&, mean radius, and the atmospheric density and its derivative at 1 bar (see details in 
Anderson & Schubert, 2007; Helled et ah, 2011; Kaspi et ah, 2013). We then use the difference in 
the values of the figure functions in the two limiting cases to limit their values. 
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Figure 9: Extended Data. The calculated flattening parameters by the model when the figure 
functions are limited by interior models, a-c A seatter plot (similar to Fig. la) of the distribution 
of solutions on the plane of the calculated rotation period Peak minus Pvoy and A J 2 , AJ 4 , AJg, 
respectively. Each blue dot represents one sub-case converged solution. In a the inner and outer 
black circles show the first and second standard deviations, respectively, d-j - the distribution of 
solutions for the figure functions K 2 , Fi, F 2 , F^, respectively. 

























Saturn 

Jupiter 

Mass (10^4 kg) 

568.36 

1,898.3 

Rotation Period (System III) 

lOhr 39m 24s 

9hr 55m 29s 

Mean Radius (km) 

58,232 

69,911 

J2 (10-6) 

16, 290.71 ±0.27 

14,696.43 ±0.21 

(10-6) 

-935.83 ±2.77 

-587.14 ±1.68 

Jg (10-6) 

86.14 ±9.64 

34.25 ± 5.22 


Table 1: Extended Data. The physical properties of Saturn and Jupiter used in the anal¬ 
ysis. Data are taken from http://ssd.jpl.nasa.gov/7gravity_fields_op, The gravitational moments 
eorrespond to a referenee equatorial radius of 60,330 km and 71,492 km for Saturn and Jupiter, 
respeetively. 











Saturn - Massive Core 

Saturn - No Core 

Fi 

0.66357 

0.67720 

F 2 

0.07501 

0.04743 

Fs 

0.043517 

0.10790 

K 2 

0.22691 

0.07535 

K 3 

-0.01039 

-0.01174 

H 3 

0.10074 

0.23489 


Jupiter - Massive Core 

Jupiter - No Core 

Fi 

0.77014 

0.76965 

F 2 

0.08512 

0.08499 

F 3 

0.14874 

0.15137 

K 2 

0.05868 

0.05861 

K 3 

-0.00705 

-0.01068 

H 3 

0.19885 

0.20570 


Table 2: Extended Data. The calculated figure functions based on interior models of Saturn 

and Jupiter. The values of the figure funetions are derived for the two limiting eases of massive 
eore and no-eore (eontinuous density profile) for Saturn and Jupiter. 



